function DerwentAWQSurfMap(varname,awqsites,data,pos,time,setasfile,ftype)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIS FUNCTION RETURNS A MAP OF AN AWQ VARIABLE IN THE DERWENT
% ESTUARY. THE MAP IS PRODUCED  BY FIRST INTERPOLATING THE SCATTERED DATA
% ONTO A GRID USING THE DELAUNAY TRIANGULATION SCHEME.
%
% INPUT: - varname  : name as in AWQ_parameters.xls  
%        - awqsites : list of awq site names
%        - data     : array of a single value per site. 
%                     Note : 1) the data must be listed as in awqsites  
%                            2) the array must be  specified in column
%                               vector format (site x value) for  correct
%                               input to the interpolation function. 
%        - pos      : position (string). Available options are 's' or 'b'
%        - time     : sampling date (matlab format)
%        - setasfile: a SETAS monthly file to provide grid informations
%        - ftype    : string, indicate type of file providing grid
%                     informations; available options are 'nc' for NetCDF 
%                     files or 'mat' if Matlab mat file
%
%
% EXAMPLE : load testmatAWQ.mat
%           DerwentAWQSurfMap('Chla_lab',get_infoAWQ('list'),...
%           testmatAWQ.Chla_lab,'b',731580,'SETAS_lon_lat.mat','mat')
%
% Author: Benedicte Pasquer, IMOS/eMII (http://imos.org.au/)
% email: benedicte.pasquer@utas.edu.au
% May 2013   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%GET SITE POSITIONS
for i =1: length(awqsites)
    
    [lo,la] = get_infoAWQ('pos',awqsites{i});
    slon(i) = lo;slat(i) = la;
     
end
%READ SETAS GRID INFORMATION 
switch ftype
    case 'nc'
% INFORMATION FROM SETAS NETCDF MONTHLY FILE
        ncid = netcdf.open(setasfile,'NC_NOWRITE');
        varid = netcdf.inqVarID(ncid,'latitude');
        lat = netcdf.getVar(ncid,varid); 
        varid = netcdf.inqVarID(ncid,'longitude');
        lon = netcdf.getVar(ncid,varid); 
        netcdf.close(ncid)
    case 'mat'
        load(setasfile)
end


%INTERPOLATION ON SETAS GRID USING DELAUNAY TRIANGULATION
%NEED TO PROVIDE ROW VECTOR OF VALUE FOR TRISCATTERED INTERPOLATION

F = TriScatteredInterp(slon',slat',data,'nearest');
zi = F(lon,lat);

%MAP
ax1 = axes('Position',[.2 .2 .65 .7 ]);
m_proj('mercator','long',[147.00 147.5],'lat',[-43.07 -42.7])
h = m_pcolor(lon,lat,zi);shading interp;
set(h,'Parent',ax1)
hold on
%DISPLAY SITE LOCATION
m_scatter(slon,slat,15,'MarkerEdgecolor','k','Linewidth',2)

%GRID DISPLAY PARAMETER
m_grid('xtick',[147 147.16667 147.3333 147.5],...
	   'ytick',[-43 -42.8333 -42.9166 -42.75],...
       'linewi',1,'linest',':',...
	   'tickdir','out','fontsize',12);
   
%USE THE DERWENT COAST DATASET INSTEAD OF BUILT IN COASTLINE DATA   
m_usercoast('Derwent.mat','patch',([1 1 1]));

%GET INFO FOR TITLE 
[longname4labl,unit] = get_longname(varname);
switch pos
    case 's'
        t1 ='surface';
    
    case 'b'
        t1 ='Bottom';  
end
% obsdate = tobs(max(~isnan(unique(tobs)))); %take last sampling day if several

titre{1} = ['Observed ',t1,' ',char(longname4labl),'  '];
titre{2} = ['in the Derwent estuary',' [',char(unit),'] on ', datestr(time,'mmm dd yyyy'), ' '];

title(titre,'Fontsize',11,'Fontweight','bold')

%CAXIS MIGHT NEED TO BE ADJUSTED ACCORDING TO DATA 
caxis([min(min(data))*0.75 max(max(data))*1.25]) 
colorbar

%CREATE A SECOND AXIS FOR FOOTNOTE
ax2 = axes('Position',[.21 .08 1 1],...
'XAxisLocation','bottom',...
'YAxisLocation','left',...
'Color','none',...
'XColor','k','YColor','k','Visible','off');
footnote{1} ='Note that this diagram is an example of potential product that could be provided.';
footnote{2}='The spatial interpolation of scatterred observations in the Derwent estuary is';
footnote{3}='produced using a Delaunay triangulation scheme. The reliability of the information';
footnote{4}='provided depends on  the main properties of the collection process, ie data quality';
footnote{5}='and spatial resolution of sampling';
text(0,0, footnote ,...
    'Color',[0 0 .7], ...
    'Fontname','arial','Fontsize',8,'Fontweight','bold',...
    'Parent',ax2)

figname =['MAP' varname pos]; 
saveas(gcf,figname,'png')

